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Abstract: We consider a hybrid array composed of scintillation and water Cherenkov detectors designed to measure the 
cosmic ray primary mass composition at energies of about 1 EeV. We have developed a simulation and reconstruction 
chain to study the theoretical performance of such an array. In this work we investigate the sensitivity of mass composition 
observables in relation to the geometry of the array. The detectors are arranged in a triangular grid with fixed 750 m 
spacing and the configuration of the scintillator detectors is optimized for mass composition sensitivity. We show that the 
performance for composition determination can be compared favorably to that of Xm^x measurements after the difference 
in duty cycles is considered. 
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1 Introduction 

The measurement of the mass-composition of ultra-high 
energy cosmic rays is one of the keys that can help us 
elucidate their origin. Up to energies around 10^^ eV, the 
Galaxy is believed to be the source of cosmic rays. Sev- 
eral acceleration mechanisms are certainly at play but it 
is widely expected that the dominant one is first order 
Fermi acceleration at the vicinity of supernova remnant 
shock waves. These Galactic accelerators should theoreti- 
cally become inefficient between 10^^ eV and lO'^ eV. The 
KASCADE experiment has measured the energy spectra 
for different mass groups in this energy range and found 
that there is a steepening of the spectra at an energy that 
increases with the cosmic ray mass [1]. As a result, the 
mass-composition becomes progressively heavy. It is also 
thought that extra-galactic sources can start to contribute 
to the total cosmic ray flux at energies above 4 x 10'^ eV. 
The onset of such an extra-galactic component would prob- 
ably produce another change in composition. The X^ax 
measurements from the HiRes-MIA experiment have been 
interpreted as a change in composition, from heavy to light, 
starting at 4 x 10'^ eV and becoming proton-dominated 
at 1.6 X lO'^eV [2, 3] while the measurements from the 
Pierre Auger observatory hint at a light or mixed composi- 
tion that becomes heavier beyond 2 X 10'^ eV [4]. 

1.1 Measuring Composition 

Roughly speaking, the techniques for inferring the mass 
composition of cosmic rays can be split in two categories, 
depending on whether they exploit the sensitivity to the 



depth of shower maximum (X^ax) or to the ratio of the 
muon and electromagnetic components of the air shower 
[5]. Direct measurements of the fluorescence emission fall 
in the first category, and so do the various measurements 
of the Cherenkov light produced by air showers. Most 
ground-based detector observables depend one way or an- 
other on the number of muons in the air shower However, 
the arrival time profile of shower particles has been used as 
an observable mostly sensitive to X^ax, in particular the so- 
called rise-time, the time it takes for the signal to rise from 
10% to 50% of the integrated signal [6]. The measurement 
of the number of muons and electrons in the air shower can 
be done directly, for example, the way it was done with the 
KASCADE detector [7]. 

The Pierre Auger Observatory is developing a series of en- 
hancements that aim at the energy range between 10'^ eV 
and 10'"^ eV [8, 9]. In particular, the objective of the 
AMIGA enhancement [8] is to measure the muon compo- 
nent of the air shower using scintillators shielded by several 
meters of soil. In the same spirit, we are considering a hy- 
brid surface array, consisting of two super-imposed ground 
arrays, a Water Cherenkov Detector (WCD) array and a 
scintillation detector array. The purpose of the scintillation 
detectors is to increase the sensitivity to the electromag- 
netic part of the air shower. 

We will consider an array of stations arranged in a triangu- 
lar grid with a separation of 750 m. Each scintillator sta- 
tion is made with 3 cm thick plastic scintillator tiles like 
the ones used for the muon detectors in the KASCADE ar- 
ray. In order to enhance the signal from gamma rays in air 
showers, we study the effect of adding a certain amount of 
lead on top of the scintillators. For conversion of around 
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80% of the high energy gamma rays one normally needs a 
shielding of about 2 radiation lengths. One radiation length 
corresponds to 0.56 cm of lead. 

1.2 Hybrid Surface Detector Performance 

In order to study such a hybrid detector, we have imple- 
mented a simulation and reconstruction chain based on the 
Pierre Auger Observatory offline framework [10]. Using 
this chain, we have processed 20,000 simulated air show- 
ers produced with the CORSIKA [11] event generator. All 
these showers were generated using QGSJET II for high 
energy hadronic interaction simulations [12]. 

The parameter space chosen for the simulations has been 
determined by the energy and angular ranges under con- 
sideration. Additionally, it has been shown recently that 
real air-showers appear to have a larger muon component 
than simulated air-showers [13]. For this reason we have 
also performed the simulations with the number of muons 
artificially increased by a factor of two, which gives the 
observed muon contribution after an energy shift of about 
28%. We have also investigated the effect of placing a cer- 
tain amount of shielding on top of each scintillator station, 
as well as the effect of changing its area. In total there are 
600 parameter combinations, with 100 simulated showers 
each, for a total of 60,000 simulated events. The parame- 
ters are shown in table 1 . 



Primaries p, Fe 

E(eV) 10'"^ 10iM0i»-25 

e 0°, 25°, 38°, 48°, 60° 

Area 3 m^, 10 m^ 

^scale 2 

Shield (rad. lengths) 0, 0.25, 0.5, 1, 2 



Table 1 : Parameters considered in the simulation of the re- 
sponse of the scintillator array. 

The simulation of the interactions of the shower parti- 
cles with the detector is done using the Geant4 package 
[14, 15]. The scintillation efficiencies used in the simula- 
tion correspond to the specifications for Bicron's BD-416 
scintillators: Polyvinyl Toluene scintillators with a nomi- 
nal light yield of about lO'* photons/Me V and a density of 
1.032 g/cm^. The resulting scintillation photons are sam- 
pled at a frequency of 100 MHz to produce one FADC trace 
per station. The signal in each station is measured in units 
of Minimum Ionizing Particle equivalent, or MIP, where a 
MIP is given by the position of the peak of the Landau dis- 
tribution for vertical muons. We then consider only stations 
with signals between 1 and 2000 MIP in order to simulate 
the dynamical range. 

The arrival direction and core position of each event are 
currently estimated using only the WCDs. The arrival di- 
rection is determined by fitting a spherical shower front 
to the signal start times of the stations in the event. The 
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Figure 2: A typical lO"* eV shower arriving at 38° showing the 
reconstructed LDF of the scintillator array. 

core position is determined by adjusting a lateral distribu- 
tion function (LDF) of the form 

Swcoir) = Sro (1) 

to the total signal in the stations in the event. The ro param- 
eter is 450 m and ^450 is the usual energy estimator for a 
WCD array with this grid spacing because it is the optimum 
distance for determining the signal [16]. The LDF for the 
scintillator array is described by the same formula and the 
optimum ro that will enhance the composition sensitivity 
still needs to be determined. 

One typical event is displayed in figures 1 and 2. An ex- 
ample of the resulting average lateral distribution functions 
over 100 showers, for both the scintillator and the WCD 
array, can be seen in figure 3. Clearly, the slope of the 
proton and iron LDFs differ and there is a cross-over point 
where the signals are equal. This cross-over point depends 
on energy and zenith angle and this dependence is differ- 
ent for the two arrays. At 10'^ eV and 0°, it is close to the 
shower axis for both arrays. As the zenith angle increases, 
the cross-over occurs at greater distances to the shower axis 
and for angles greater than 25°, which correspond to most 
of the aperture, it is located at more than 600 meters from 
the axis. At larger distances, the differences between pro- 
ton and iron signals are similar in the scintillators and the 
WCDs since they both detect mostly muons in this radial 
range. Therefore, the largest difference in signal between 
proton and iron is found close to the axis. 

For the reason mentioned above, we use the signal at 400 
meters from the shower axis as the measure for the scin- 
tillator array signal. For each event we can then correlate 
the scintillator array signal (S'sci) with the WCD signal at 
450 m (S'wcd)- A typical contour plot of these quantities 
is shown in figure 4. One can model the dependence of the 
scintillator and WCD signals as a function of the energy 
and the mass number using 

logio S, = a, + b, logio E + c, logio ^ (2) 
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Figure 1 : Comparison of WCD and scintillator traces. These two sets of traces correspond to a pair of scintillator station with no 
photon converter (left) and a WCD station (right). They are located at 413 m from the axis of a lO' eV shower arriving at 38°. Note 
that, while the WCD station already has a significant contribution from muons, this is very small in the scintillator, which is dominated 
by the signal from electrons. 
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Figure 3: Comparison of the lateral distributions of 1 x lO'^ eV 
proton and iron induced air showers arriving at 38 degree using 
a combined scintillatorAVCD array. The area of each station is 
10 m~. The muon component is artificially scaled by a factor of 
two. 

This set of equations can then be inverted to find logj^Q E 
and logiQ A. The result is a mass estimator and an unbiased 
energy estimator. 

There are various ways of measuring the discrimination 
power of a given statistic. One of them is the figure of 
merit: 

f=i^fi=M. (3) 

The resulting average discrimination power, after correct- 
ing by the decrease in aperture with increasing zenith angle, 
is displayed in figure 5. We have used other multivariate 
methods as well as alternative measures of discrimination 
power, such as the efficiency at given purity levels. The 
results are qualitatively the same. 

As an exercise to get an idea of the magnitude of the result- 
ing uncertainties in the determination of the composition, 
one can pose the problem of determining the fraction of 
protons in a sample of proton and iron events. 

Given a sample of equal number proton and iron showers, 
we can determine a cut value for the mass estimator such 



Figure 4: 68% contours of thelogj,)(5'sci) -log^g(S'wcD) scatter 
plot. This corresponds to a zenith angle of 25°, 10 m^ stations 
with no lead photon converter and a scaling of the muon compo- 
nent by a factor of two. 

that the error rates of the first and second kind are equal (e). 
Let's now consider a sample of N events where a fraction 
/ of them are protons. The observed number of /?roton-like 
{Sobs, the number of events passing the cut) and /ron-like 
(bobs) events depend on the number of actual signal and 
background events through the following equation: 

( bl ) = ( e 1-e) { (l-m) 

and this equation can be inverted to estimate the proton 
fraction. 

In order to judge the quality of the method, a common strat- 
egy is to split the sample into a training sample, from which 
we get the values of the coefficients in equations 2 and the 
value of the cut. We can also use a k-fold cross validation 
method. This consists in splitting the sample in sets of k 
events, taking half of the sets to form the training sample 
and use the rest as a test sample. This process is then re- 
peated with different combinations of sets. The same pro- 
cedure is repeated using the distributions of v^Cmax for pro- 
ton and iron but in this case we limit the size of the test 
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Figure 5; Average figure of merit as a function of tiie thickness 
of tiie lead siiield. Tiie different graphs correspond to the different 
detector areas and muon component scaling factor considered. 




Figure 6: Cross validation of the proton fraction determination. 
This corresponds to lO'** eV showers arriving at 38° and muon 
component scaled by a factor of two. The stations have an area of 
3 m^ and 0.25 radiation length photon converters. The test sample 
for the Xmax technique is 30% of the sample for the combined 
array. 



sample to 30% in order to take the combined effect of the 
aperture and limited duty cycle of a fluorescence detector 
into account. The result can be seen in figure 6. 

From this simple exercise we have learned that both meth- 
ods give similar uncertainties in the determination of the 
proton fraction. In other words, that the combined scin- 
tillator/WCD array provides a composition estimate that is 
competitive with that from fluorescence detectors after one 
considers the fluorescence detectors' limited duty-cycle. 
We are aware that the precise factor to penalize the X^ax 
measurement will depend on the relative apertures of the 
fluorescence detector and the surface detector and that there 
will be systematic errors deriving from the different strate- 
gies in training the classification and regression methods. 
This number wiU still give the right order of magnitude. 



2 Summary 

We have performed detailed simulations of the response of 
a combined scintillator/WCD array to air showers produced 
by cosmic ray with primary energies around 1 x 10'^ eV. 
We have considered different configurations for the indi- 
vidual scintillator stations in order to optimize the compo- 
sition sensitivity of the combined array. 

We have found that an array of lOm^ scintillator stations, 
together with the WCD aiTay, can provide an estimate of 
the average mass composition that is competitive with X,nax 
measurements after differences in duty-cycle are consid- 
ered and that the addition of photon converters on top of 
each station does not add to the separation power. 

In this work we only considered a grid spacing of 750 m 
but we have noted that the highest separation between pro- 
ton and iron occurs at small distances to the core. Better 
knowledge of the impact of smaller spacing and of the pre- 
cise optimum distance at which to measure the scintillator 
signal requires further study. 
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